%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IMPORTS HERE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass{beamer}
\usetheme{Warsaw}
\usepackage[utf8]{inputenc}
\usepackage{default}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{tikz}
\usepackage{graphicx}
\usepackage{mdframed}
\usepackage{color}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SETUP SECTION
% Declare any macros here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\mode<presentation>{
  \usetheme{Warsaw}
  \setbeamercovered{transparent}
}
\AtBeginSection[]{
  \begin{frame}<beamer>{Outline}
    \tableofcontents[currentsection,currentsubsection]
  \end{frame}
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% AUTHOR AND TITLE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{Notes on imaging}
\author{Benjamin Hugo}
\date{March 2014}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% START OF PRESENTATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
% \begin{frame}
%  \maketitle
% \end{frame}
% \section{Overview of the process (sunny day)}
% \begin{frame}{Sampling on the celestrial sphere}
%  \begin{center}
%   \includegraphics[height=0.8\textheight]{Images/cel_sphere.png}
%  \end{center}
% \end{frame}
% 
% \begin{frame}{Sampling at phase center}
%  \begin{center}
%   \includegraphics[height=0.9\textheight]{Images/sampling_coherence.png}
%  \end{center}
% \end{frame}
% 
% \begin{frame}{The lmn space}
%  \begin{center}
%   \includegraphics[width=0.7\textwidth]{Images/lmn_to_uvw.png}
%  \end{center}
% \end{frame}

\section{Problems in imaging, deconvolution and calibration}

\begin{frame}{Resolution}
 \begin{itemize}
  \item<1-> improves with larger dishes
  \item<2-> improves with longer baselines
  \item<3-> decreases for large $\lambda$ (low frequency)
 \end{itemize}
\end{frame}

\begin{frame}{from uvw to lmn and back - the measurement equation (TACA)}
  \footnotesize
  \begin{equation*}
    V_{pq}(\sigma,t,\upsilon) = G\int \mathbf{D_p}K_p \left[ \begin{array}{cc} e_{XX} & e_{XY} \\ e_{YX} & e_{YY} \end{array} \right] K_q^H\mathbf{D_q}^H d\Omega G^H
  \end{equation*}
    where: 
  \begin{itemize}
    \item $D_p,K_p,B,K_q^H,D_q^H$ are functions of $(\sigma,t,\upsilon)$. $G,G^H$ are source independent.
    \item $e_{XX}, e_{XY}, e_{YX}, e_{YY} \in \mathbb{C}$ the correlations between 2 polarizations
    \item $K_p K_q = e^{\frac{-2\pi i}{\lambda} \vec{b}(\vec{s}(\upsilon) - \vec{s_0}(\upsilon))} I$, where $I$ is the 2x2 identity matrix and $\vec{b}(\upsilon)(\vec{s}(\upsilon) - \vec{s_0}(\upsilon))=u_{pq}l+v_{pq}m+w_{pq}(\sqrt{1-l^2-m^2}-1)$ for each $\upsilon$
    \item \textbf{If $\sqrt{1-l^2-m^2} \approx 1$ the integral is the fourier transform between the \underline{dirty} image and the coherence function.}
  \end{itemize}
\end{frame}

\begin{frame}{Countermeasures!}
\Tiny
 \begin{itemize}
  \item \textbf{Directional dependent effects}\only<2>{ include ionospheric distortions and Faraday rotations. Calibration required to obtain the inverse matrices for each $D(\sigma,t,\upsilon)$}
  \item \textbf{Incomplete sampling in the fourier domain.}\only<3>{ Causes aliasing in the image (\textit{PSF}). Countermeasure: \textit{Deconvolution}.}
  \item \textbf{Time \& Frequency smearing (radial).}\only<4>{ 
  \begin{center}
      \includegraphics[width=0.6\textwidth]{Images/radial_smearing.png}
  \end{center}
  Fringe modulation far from the delay center. Countermeasure: limit observation bandwidth and average coherence function measurements must be over short intervals: $<e_pe_q^*>$ 
  }
  \item \textbf{An single reflector telescope is sensitive to a limited field of view $\approx\frac{\lambda}{D}$. An interferometer is sensitive to $\approx\frac{\lambda}{B}$} \only<5>{ Countermeasure: mosaicking technique required for large field of view.}
  \item \textbf{Can only use the measurement equation if $\sqrt{1-l^2-m^2} \approx 1$ (only over a small FOV).} \only<6> { Countermeasure: faceting }
  \item \textbf{PSF stays relatively constant only over a small FOV.} \only<7> { Countermeasure: faceting + CLEAN}
  \item \textbf{The problem of co-planar sampling: Fresnel distortion $\propto\Delta w$ .} \only<8>{ 
    \begin{center}
      \includegraphics[width=0.3\textwidth]{Images/fresnel_distortion.png}
    \end{center}
    Countermeasure: faceting or w-projection or snapshots algorithm
  }
 \end{itemize}
\end{frame}
\section{Simple pipeline}
\begin{center}
\begin{frame}{Imaging and deconvolving}
 \begin{tikzpicture}[font=\tiny]
    \draw (0,0) rectangle (2,1) node at (1,0.5) {Antennae};
    \draw[->] (2,0.5) -- (3,0.5);
    \draw (3,0) rectangle (6,1) node at (4.5,0.5) {Interferometer};
    \draw[->] (6,0.5) -- (7,0.5);
    \draw (7,0) rectangle (9,1) node at (8,0.5) {$V - V_{model}$};
    \draw[->] (8,1) -- (8,2);
    \draw (7,2) rectangle (9,3) node at (8,2.5) {Gridding};
    \draw[->] (8,3) -- (8,4) node[pos=.5, right]{Dirty $V_i^s$};
    \draw (7,4) rectangle (9,5) node at (8,4.5) {IFFT};
    \draw[->] (7,4.5) -- (6,4.5) node[pos=.5, above]{Dirty $I_i^s$};
    \draw (3,4.5) -- (4.5,4) -- (6,4.5) -- (4.5,5) -- cycle node at (4.5,4.5) {Converged};
    \draw[->] (3,4.5) -- (2,4.5);
    \draw (0,4) rectangle (2,5) node at (1,4.5) {Deconv. \& Update};
    \draw[->] (1,4) -- (1,3) node[pos=.5, right]{$I_{model}^s$};
    \draw (0,2) rectangle (2,3) node at (1,2.5) {FFT};
    \draw[->] (2,2.5) -- (3,2.5) node[pos=.5, above]{$V_{model}^s$};
    \draw (3,2) rectangle (6,3) node at (4.5,2.5) {Degridding};
    \draw[->] (6,2) -- (7,1) node[pos=.5, sloped, above]{$V_{model}$};
    \draw[->] (4.5,5) -- (4.5,6) node[pos=.5, right]{CLEAN $I_m^s$};
    \draw[style=thick] (4.4,6.1) -- (4.6,6.1);
    \draw[style=thick] (4.3,6.2) -- (4.7,6.2);
 \end{tikzpicture}
\end{frame}
\end{center}
\begin{frame}{Gridding}
 \tiny
 Problem:
 \begin{itemize}
  \item<1-> UV plane is continuous. Taking an FFT / IFFT requires uniformly spaced samples.
  \item<2-> \textbf{Sampling it will likely create aliasing artifacts!}
  \item<3-> Need to use a \textbf{linearly separable} anti-aliasing FIR. (convolution with infinite \textit{sinc} function achieves optimal gridding --- obv. not possible). Prolate 
    Spheriodal Wave Functions used instead (TBA).
  \item<4-> Need to grid XX,XY,YX and YY correlations \textbf{separately}
 \end{itemize}
\end{frame}
\begin{frame}{Gridding illustration (3x3 FIR). Steps 1 - 2} 
 \begin{tikzpicture}[font=\tiny]
  \filldraw[fill=red!20] (4.50,1.75) rectangle (5.25,2.5);
  \filldraw[fill=blue!20] (3.75,2.75) rectangle (4.50,3.50);
  \filldraw[fill=green!20] (5.50,4.75) rectangle (6.25,5.5);
  \only<2>{
  \draw[step=0.25,gray,thin] (3,1.5) grid (7,5.5);
  }
  \draw[red] node at (0,0.25) {data at continuous coords (per baseline, per channel):};
  \draw node at (6.5,0.25) {\dots};
  \draw[step=1,black] (3,0) grid (6,0.5);
  \draw node at (3.5,0.25) {$(u_0,v_0,\mathbb{C})$};
  \draw node at (4.5,0.25) {$(u_1,v_1,\mathbb{C})$};
  \draw node at (5.5,0.25) {$(u_2,v_2,\mathbb{C})$};
  \draw[->] (3.5,0.5) -- (5-0.125,2.25-0.125);
  \draw[->] (4.5,0.5) -- (4.25-0.125,3.25-0.125);
  \draw[->] (5.5,0.5) -- (6-0.125,5.25-0.125);
  \draw[->,gray] (3,1.5) -- (8,1.5) node[below,pos=0.5]{u};
  \draw[->,gray] (3,1.5) -- (3,6.5) node[left,pos=0.5]{v};
 \end{tikzpicture}
\end{frame}
\begin{frame}{Gridding Algorithm}
    \tiny
    \begin{enumerate}
     \item The V(u,v) data is sampled ($V(u,v)\cdot PSF(u,v)$) according to some irregular 2D sampling function:\\
	   \begin{equation*}
	    PSF_\epsilon(u,v) = \sum_{j=1}^{P}{\delta(u-u_j,v-v_j)}\cdot\epsilon(u,v)
	   \end{equation*}
	   $\epsilon(u,v)$ may be a weight (as with natural weighting) or a function of how many sampling points lie in each grid cell (uniform weighting).
     \item The data is then convolved with C and sampled with the comb function $III(u,v)$:\\
	  \begin{equation*}
	  V^s(u,v) = ([V(u,v)\cdot PSF(u,v)]*C(u,v)) \cdot III(u,v)
	  \end{equation*}
	  Where the comb function is defined as:\\
	  \begin{equation*}
	    III(u,v) = \sum_{i=-\infty}^{\infty}\sum_{j=-\infty}^{\infty}\delta(u-i\Delta u,v-j\Delta v)
	  \end{equation*}
     \item Take $\mathcal{F}^{-1}\{V^s(u,v)\}$:
	  \begin{equation*}
	    I^s_{dirty} = ([I(l,m)*\mathcal{F}^{-1}\{PSF(u,v)\}]\cdot\mathcal{F}^{-1}\{C(u,v)\}) * \mathcal{F}^{-1}\{III(u,v)\}
	  \end{equation*}
     \item To avoid attenuation at the edges introduced by C(u,v) with non-zero rolloff we divide through by $\mathcal{F}^{-1}\{C(u,v)\}$ and to discard replicas (at intervals 
	   $((\Delta u)^{-1},(\Delta v)^{-1})$) created by convolution with the $W\times W$ dimension comb function we multiply by the boxcar function:
           \begin{equation*}
	    \sqcap (l,m) = \left\{
				  \begin{array}{ll} 
				  1 & if |l| < 0.5W \wedge |m| < 0.5W \\ 
				  0 & otherwise 
				  \end{array}
			  \right.
           \end{equation*}
    \end{enumerate}
\end{frame}

\begin{frame}{Gridding: some notes} 
\tiny
\begin{itemize}
 \item<1> The PSF must be compact: if there are too few visibilities the PSF becomes excessively large and cannot be deconvolved.
 \item<2> The PSWF functions attempt to minimize aliasing energy introduced by rolloff in a typical bandpass filter. These functions and their transforms are computationally intensive, but 
    the Kaiser-Bessel function is a good approximation to these PSWFs:\only<2>{
    \begin{equation*}
      C(u) = \frac{1}{W}B_0[\beta\sqrt{1 - \left(\frac{2u}{W}\right)^2}],\beta\text{ minimizes }J = \frac{\int_{|u|>0.5W}{\left| \mathcal{F}^{-1}\{C(u)\}\cdot\left[ \frac{\sqcap(u)}{\mathcal{F}^{-1}\{C(u)\}} * \mathcal{F}^{-1}{III(u)} \right] \right|^2 du}}{\int_{-\infty}^{\infty}\left| \mathcal{F}^{-1}\{C(u)\}\cdot\left[ \frac{\sqcap(u)}{\mathcal{F}^{-1}\{C(u)\}} * \mathcal{F}^{-1}{III(u)} \right] \right|^2 du}
    \end{equation*}
    Where $B_0$ is the Bessel function of order 0:
    \begin{equation*}
      B_0 = \sum_{n=0}^{\infty}\frac{(-1)^n x^{2n}}{2^{2n}(n!)^2}
    \end{equation*}
    The inverse transform of the Kaiser-Bessel is:
    \begin{equation*}
      \mathcal{F}^{-1}\{C(u)\} = \frac{\sin{\sqrt{\pi^2 W^2 u^2 - \beta^2}}}{\sqrt{\pi^2 W^2 u^2 - \beta^2}}
    \end{equation*}
    }
 \item<3> The Keiser-Bessel gives a good approximation to a PSWF (limiting aliasing energy).\only<3>{
    \begin{center}
      \includegraphics[width=0.45\textwidth]{Images/keiser-bessel.png}
    \end{center}
 }
 \item<4> The grid sampling rate must comply with Nyquist.\only<4>{
    \begin{equation*}
      \Delta u = \frac{1}{2|u_{max}|} \wedge \Delta v = \frac{1}{2|v_{max}|}
    \end{equation*}
 }
 \item<5> The convolution function is oversampled by a factor of 10-100 and must have support preferably over the entire uv space.\only<4>{
    \begin{equation*}
      \beta_u = \frac{\Delta u}{\Delta u_{conv}} \wedge \beta_v = \frac{\Delta v}{\Delta v_{conv}}
    \end{equation*}
 }
\end{itemize}
\end{frame}

% \begin{frame}{Degridding}
%  \begin{itemize}
%   \item<1> Reverse the correcting steps in step 4 of the gridding algorithm
%   \item<2> Over all the visibilities in time (per baseline and channel) grab the contributions from all cells within the support area of $C(u,v)$ centered
%   at the real valued coordinates $(u,v)$. Weigh them with $C(u,v)$ and accumulate all the coordinates from the grid in the support area.
%   \item<3> Normalize the value of the coordinate and store.
%  \end{itemize}
% \end{frame}

\begin{frame}{W-projection}
 \tiny
 \begin{itemize}
  \item<1> Corrects the Fresnel distortion and can be done in uvw or image space.\only<1>{
  It can be shown that any plane $V_w$ can be projected to $V_0$ by convolving with a w-dependent function $G^{-1}_w(u,v)$:
    \begin{equation*}
     V_0(u,v) = V_w(u,v) * G_w^{-1}(u,v)
    \end{equation*}
    Where the function is invertible in image space and is defined as:
    \begin{equation*}
      \mathcal{F}^{-1}\{G(u,v)\} = e^{-2\pi iw(\sqrt{1-l^2-m^2}-1)}
    \end{equation*}
  }
  \item<2> Can be integrated into the gridding algorithm (requires multiple oversampled convolution functions over various regions per value of w - think 
    of this as a 3D array). This likely requires a lot less memory than doing w-stacking in image space;
  \item<3> w-stacking is done by gridding multiple planes (where the w coordinates in uvw space snap to the closest plane). We then multiply each plane 
    with $\mathcal{F}^{-1}\{G^{-1}(u,v)\}$ and merge them with $I_0^s$.
  \item<4> w-stacking is claimed to be computationally faster than convolutional w-projection. w-projection using snapshots has been shown to be faster 
    than convolutional w-projection.
  \item<5> w-stacking is an order of magnitude faster than faceting.
 \end{itemize}
\end{frame}

\begin{frame}{The Polyhedron Method (in uvw space)}
\tiny
 \begin{itemize}
  \item<1> Break the celestial sphere up into many small rectangular planes where $n \approx 1$ and the Fresnel distortion can be ignored. All directional effects \& the PSF are 
  assumed constant and can be treated after degridding and transform.\only<1>{To eliminate phase error we need:
  \begin{equation*}
    N_{\text{facets per dim}} = \frac{2\lambda B_{max}}{fD^2}
  \end{equation*}
  Where f = maximum number of pixels from sphere ($\sim 0.2$ for HDR imaging).
  }
  \item<2> Algorithm involves a phase shift and a rotation of the basis uvw vectors\only<2>{
  \begin{enumerate}
    \tiny
    \item Begin by phase shifting all the visibilities to the new phase center on the unit sphere:
    \begin{equation*}
     (\forall u,v,w) V^{PR}(u,v,w) = V(u,v,w)e^{2\pi i[u(l_0-l)+v(m_0-m)+w(n_0-n)]}
    \end{equation*}
    {\color{red} Assumption: initial $(l_0,m_0,n_0)$ is the intersection of $\vec{w}$ with the unit sphere, extending $\frac{N_{\text{facets per dim}}}{2}$ in $\pm l,m$ directions. 
    Compute $n_0-n$ from unit sphere equation.}
    \item Rotate the basis vectors uvw so that w points to the new phase center on the unit sphere $(r,\Delta\theta,\Delta\phi)$ in spherical coordinates). Thereafter do a rotation on
      $\vec{u},\vec{v},\vec{w}$. This step is critical to avoid a decrease in undistorted FOV.
      \begin{center}
      \includegraphics[width=0.3\textwidth]{Images/no_uvw_transform.png}
      \end{center}
    \item Grid each facet separately.
    \item {\color{red} Papers only briefly mentions problem of overlap in deconvolution stages (esp. degridding)}
    \item {\color{red} Image space faceting is ``depricated'' as Cornwell puts it (apparently leads to discontinuities in image?)}
  \end{enumerate}
  }
 \end{itemize}
\end{frame}
% \begin{frame}
% \color{red}
%  Where to from here (TODO):
%  \begin{itemize}
%   \item Read Synthesis Imaging Taylor, Carilli \& Perley - including especially calibration and CLEAN.
%   \item Gather more info on snapshots.
%   \item Need to look into xyz $\Longleftrightarrow$ uvw $\Longleftrightarrow$ ra/dec coordinate transforms.
%   \item Read Oleg's papers on RIME-based calibration ($3^{rd}$ gen).
%  \end{itemize}
% 
% \end{frame}


\end{document}

